
set.seed(12345)
########## Model 1 ##############

z.out1 <- zelig(non_comply_dummy ~ congress_required,model ="logit", data=data2)
print(summary(z.out1), digits = 2)

x.high1 <- setx(z.out1, congress_required=1)
x.low1 <- setx(z.out1, congress_required=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)


0.42, 0.18, 0.014,   0.7,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA



########## Model 2 ##############

z.out1 <- zelig(non_comply_dummy ~ congress_required+divided_government+cow_type1_dummy+usa_exports_ln+gdppc_ln+population_ln+polity2+trade_remedy_dummy,model ="logit", data=data2)
print(summary(z.out1), digits = 2)


x.high1 <- setx(z.out1, congress_required=1)
x.low1 <- setx(z.out1, congress_required=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

x.high1 <- setx(z.out1, divided_government=1)
x.low1 <- setx(z.out1, divided_government=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

x.high1 <- setx(z.out1, cow_type1_dummy=1)
x.low1 <- setx(z.out1, cow_type1_dummy=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$usa_exports_ln)
x.high1 <- setx(z.out1, usa_exports_ln=13.30)
x.low1 <- setx(z.out1, usa_exports_ln=4.847)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$gdppc_ln)
x.high1 <- setx(z.out1, gdppc_ln=12.020)
x.low1 <- setx(z.out1, gdppc_ln=6.009)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$population_ln)
x.high1 <- setx(z.out1, population_ln=21.57)
x.low1 <- setx(z.out1, population_ln=11.31)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$polity2)
x.high1 <- setx(z.out1, polity2=10)
x.low1 <- setx(z.out1, polity2=-7)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

x.high1 <- setx(z.out1, trade_remedy_dummy=1)
x.low1 <- setx(z.out1, trade_remedy_dummy=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)


0.58, 0.24, 0.045,  0.91,
-0.37, 0.2, -0.73, 0.00065,
0.44, 0.31, -0.57,  0.83,
 0.43, 0.48, -0.82,  0.98,
 -0.33, 0.52, -0.99,  0.82,
 -0.35, 0.55, -0.99,   0.7,
 -0.73, 0.37, -0.99,  0.35,
 0.2, 0.2, -0.21,  0.59,
 NA, NA,NA, NA



########## Model 3 ##############

z.out1 <- zelig(non_comply_dummy ~ congress_required+divided_government+cow_type1_dummy+usa_exports_ln+gdppc_ln+population_ln+polity2+trade_remedy_dummy+contribution_2010_ln,model ="logit", data=data2)
print(summary(z.out1), digits = 2)


x.high1 <- setx(z.out1, congress_required=1)
x.low1 <- setx(z.out1, congress_required=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

x.high1 <- setx(z.out1, divided_government=1)
x.low1 <- setx(z.out1, divided_government=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

x.high1 <- setx(z.out1, cow_type1_dummy=1)
x.low1 <- setx(z.out1, cow_type1_dummy=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$usa_exports_ln)
x.high1 <- setx(z.out1, usa_exports_ln=13.30)
x.low1 <- setx(z.out1, usa_exports_ln=4.847)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$gdppc_ln)
x.high1 <- setx(z.out1, gdppc_ln=12.020)
x.low1 <- setx(z.out1, gdppc_ln=6.009)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$population_ln)
x.high1 <- setx(z.out1, population_ln=21.57)
x.low1 <- setx(z.out1, population_ln=11.31)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$polity2)
x.high1 <- setx(z.out1, polity2=10)
x.low1 <- setx(z.out1, polity2=-7)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

x.high1 <- setx(z.out1, trade_remedy_dummy=1)
x.low1 <- setx(z.out1, trade_remedy_dummy=0)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)

summary(data2$contribution_2010_ln)
x.high1 <- setx(z.out1, contribution_2010_ln=19.86)
x.low1 <- setx(z.out1, contribution_2010_ln=17.62)
s.out1 <- sim(z.out1, x1 = x.high1, x = x.low1) 
print(summary(s.out1), digits = 2)


0.63, 0.26, 0.01,  0.94,
-0.44, 0.22, -0.82, -0.018,
0.42, 0.31, -0.54,  0.81,
0.35, 0.52, -0.9,  0.97,
-0.26, 0.52, -0.98,  0.85,
-0.25, 0.56, -0.99,  0.82,
-0.75, 0.35,   -1,  0.32,
0.26, 0.22, -0.24,  0.67,
-0.2, 0.35, -0.84,  0.53


##################################
######### Graphing ################
##################################


m <- 3 ######## Number of Tables
n <-9  ######## Number of Variables


par(mfrow=c(1,m), family = "serif",oma = c(2,2,4,2),mar = c(4,.5,4,.5))



######### Model 1 ################

vnames <- c("Congress Required",NA,NA,NA,NA,NA,NA,NA,NA)

vals <- c(
0.42, 0.18, 0.014,   0.7,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA,
NA,	NA,	NA,	NA
)

plot(NULL, col = "grey10",xlim = c(-1,1), ylim = c(.7, (n+.4)),  axes = F, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i")

polygon(c(-1.015, 1.015, 1.015, -1.015), c(.7,.7,(n+.4),(n+.4)), border = "grey", col = "grey", xpd=NA)
#text(-1, .7, "-1", pos = 1, xpd = T, cex = .8)
text(-.5, .7, "-0.5", pos = 1, xpd = T, cex = .8)
text(0, .7, "0", pos = 1, xpd = T, cex = .8)
text(.5, .7, "0.5", pos = 1, xpd = T, cex = .8)
#text(1, .7, "1", pos = 1, xpd = T, cex = .8)
abline(v = .5, lty = 3, col = "white")
abline(v = -.5, lty = 3, col = "white")
abline(v = 0, col = "white", lwd = 1)

fd <- matrix(nrow=n, ncol =4, vals, byrow=T)
for(i in 1:n){
	if(!is.na(fd[i,3])&fd[i,3] <0 &fd[i,4] > 0){
		mypch <-21
	} else { mypch <-19
		} 
					if(!is.na(fd[i,3])&fd[i,3] <0 &fd[i,4] > 0){
		myline <-2
	} else { myline <-1
		}
		lines(c(fd[i,3], fd[i,4]), c(((n+1)-i),((n+1)-i)), lwd=2,lty=myline)
		points(fd[i,1],((n+1)-i),pch = mypch, cex =.5, bg="white")
		text(fd[i,1],((n+1)-i),vnames[i],xpd=T,cex=.8,pos=3,font=2)
}

mtext(side = 1, expression(paste(Delta," in predicted")), line = 1.2, cex = .8)
mtext(side = 1, "probabilities", line = 1.9, cex = .8)
mtext(side = 3, "Model 1\nCompliance\n(n = 37)", line = 1, cex = .9)



######### Model 2 ################

vnames <- c("Congress Required", "Divided Government", "Formal Alliance", "USA Exports","GDP Per Capita", "Population", "Polity Score", "Trade Remedy",NA)
vals <- c(
0.58, 0.24, 0.045,  0.91,
-0.37, 0.2, -0.73, 0.00065,
0.44, 0.31, -0.57,  0.83,
 0.43, 0.48, -0.82,  0.98,
 -0.33, 0.52, -0.99,  0.82,
 -0.35, 0.55, -0.99,   0.7,
 -0.73, 0.37, -0.99,  0.35,
 0.2, 0.2, -0.21,  0.59,
 NA, NA,NA, NA
)

plot(NULL, col = "grey10",xlim = c(-1,1), ylim = c(.7, (n+.4)),  axes = F, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i")

polygon(c(-1.015, 1.015, 1.015, -1.015), c(.7,.7,(n+.4),(n+.4)), border = "grey", col = "grey", xpd=NA)
#text(-1, .7, "-1", pos = 1, xpd = T, cex = .8)
text(-.5, .7, "-0.5", pos = 1, xpd = T, cex = .8)
text(0, .7, "0", pos = 1, xpd = T, cex = .8)
text(.5, .7, "0.5", pos = 1, xpd = T, cex = .8)
#text(1, .7, "1", pos = 1, xpd = T, cex = .8)
abline(v = .5, lty = 3, col = "white")
abline(v = -.5, lty = 3, col = "white")
abline(v = 0, col = "white", lwd = 1)


fd <- matrix(nrow=n, ncol =4, vals, byrow=T)
for(i in 1:n){
	if(!is.na(fd[i,3])&fd[i,3] <0 &fd[i,4] > 0){
		mypch <-21
	} else { mypch <-19
		} 
					if(!is.na(fd[i,3])&fd[i,3] <0 &fd[i,4] > 0){
		myline <-2
	} else { myline <-1
		}
		lines(c(fd[i,3], fd[i,4]), c(((n+1)-i),((n+1)-i)), lwd=2,lty=myline)
		points(fd[i,1],((n+1)-i),pch = mypch, cex =.5, bg="white")
		text(fd[i,1],((n+1)-i),vnames[i],xpd=T,cex=.8,pos=3,font=2)
}

mtext(side = 1, expression(paste(Delta," in predicted")), line = 1.2, cex = .8)
mtext(side = 1, "probabilities", line = 1.9, cex = .8)
mtext(side = 3, "Model 2\nCompliance\n(n = 37)", line = 1, cex = .9)


######### Model 3 ################

vnames <- c("Congress Required", "Divided Government", "Formal Alliance", "USA Exports","GDP Per Capita", "Population", "Polity Score", "Trade Remedy","Contributions")


vals <- c(
0.63, 0.26, 0.01,  0.94,
-0.44, 0.22, -0.82, -0.018,
0.42, 0.31, -0.54,  0.81,
0.35, 0.52, -0.9,  0.97,
-0.26, 0.52, -0.98,  0.85,
-0.25, 0.56, -0.99,  0.82,
-0.75, 0.35,   -1,  0.32,
0.26, 0.22, -0.24,  0.67,
-0.2, 0.35, -0.84,  0.53
)

plot(NULL, col = "grey10",xlim = c(-1,1), ylim = c(.7, (n+.4)),  axes = F, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i")

polygon(c(-1.015, 1.015, 1.015, -1.015), c(.7,.7,(n+.4),(n+.4)), border = "grey", col = "grey", xpd=NA)
#text(-1, .7, "-1", pos = 1, xpd = T, cex = .8)
text(-.5, .7, "-0.5", pos = 1, xpd = T, cex = .8)
text(0, .7, "0", pos = 1, xpd = T, cex = .8)
text(.5, .7, "0.5", pos = 1, xpd = T, cex = .8)
#text(1, .7, "1", pos = 1, xpd = T, cex = .8)
abline(v = .5, lty = 3, col = "white")
abline(v = -.5, lty = 3, col = "white")
abline(v = 0, col = "white", lwd = 1)


fd <- matrix(nrow=n, ncol =4, vals, byrow=T)
for(i in 1:n){
	if(!is.na(fd[i,3])&fd[i,3] <0 &fd[i,4] > 0){
		mypch <-21
	} else { mypch <-19
		} 
					if(!is.na(fd[i,3])&fd[i,3] <0 &fd[i,4] > 0){
		myline <-2
	} else { myline <-1
		}
		lines(c(fd[i,3], fd[i,4]), c(((n+1)-i),((n+1)-i)), lwd=2,lty=myline)
		points(fd[i,1],((n+1)-i),pch = mypch, cex =.5, bg="white")
		text(fd[i,1],((n+1)-i),vnames[i],xpd=T,cex=.8,pos=3,font=2)
}


mtext(side = 1, expression(paste(Delta," in predicted")), line = 1.2, cex = .8)
mtext(side = 1, "probabilities", line = 1.9, cex = .8)
mtext(side = 3, "Model 3\nCompliance\n(n = 37)", line = 1, cex = .9)




